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The one-dimensional Holstein-Hubbard model with two electrons of opposite spin is studied using 
an extension of a recently developed quantum Monte Carlo method, and a very simple yet rewarding 
variational approach, both based on a canonically transformed Hamiltonian. The quantum Monte 
Carlo method yields very accurate results in the regime of small but finite phonon frequencies, char- 
acteristic of many strongly correlated materials such as, e.g., the cuprates and the manganites. The 
influence of electron-electron repulsion, phonon frequency and temperature on the bipolaron state 
is investigated. Thermal dissociation of the intersite bipolaron is observed at high temperatures, 
and its relation to an existing theory of the manganites is discussed. 
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I. INTRODUCTION 

In recent years, the formation and properties of bipo- 
larons, consisting of two electrons forming a pair in 
real space, have received considerable interest due to 
their potential role, e.g., in high-temperature supercon- 
ductivity. Theories based on bipolaron formation have 
been proposed for the superconducting transition in the 
cupratesji and the metal-insulator transition and colos- 
sal magnetoresistance in the manganites. ^■■^ Despite some 
fundamental problems, ^'^'^ they are still issue of ongoing 
discussion. 

Many interesting materials fall into the adiabatic 
regime of small but finite phonon frequencies and inter- 
mediate to strong electron-phonon coupling. For such 
parameters, analytical approaches based on, e.g., pertur- 
bation theory, do not give reliable results. In contrast, 
computational methods represent a very powerful instru- 
ment to obtain exact, unbiased information, and a lot of 
numerical work has recently been devoted to an under- 
standing of the Holstein and the Holstein-Hubbard (HH) 
model. 

In this paper, we present a simple but surprisingly ac- 
curate variational approach to the HH bipolaron. More 
importantly, we extend a recently developed quantum 
Monte Carlo (QMC) method^ to the case of two electrons 
of opposite spin. The resulting algorithm is used to study 
bipolaron formation in the one-dimensional HH model, 
focusing on the adiabatic regime. While the ground- 
state properties of the HH bipolaron are rather well un- 
derstood, here we exploit the capability of the QMC ap- 
proach to also study finite temperatures. We find that, in 
particular, the weakly bound intersite bipolaron is sus- 
ceptible to thermal dissociation. Furthermore, in con- 
trast to previous studies, we are able to consider a very 
large range of the electron-phonon and electron-electron 
interaction. 

The outline of this work is as follows. In Sec.|n]we dis- 
cuss the HH model with two electrons, while in Sec. IIIII 
we present an extended Lang-Firsov transformation with 
nonlocal lattice displacements. Section IIVI features the 



extension of the QMC method to the bipolaron problem, 
and Sec. rVI covers the variational approach. Results are 
presented in Sec. IVII and Sec. IVIII contains our conclu- 
sions. 

II. THE HOLSTEIN-HUBBARD MODEL 

The HH model is defined in terms of dimensionless 
phonon by the Hamiltonian 
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where K describes the hopping of electrons, P corre- 
sponds to the sum of the kinetic {Pp) and elastic {Px) 
energy of the phonons, and I^p, Ice denote the electron- 
phonon (el-ph) and electron-electron (cl-el) interaction 
terms, respectively. Here (cio-) creates (annihilates) 
an electron of spin a at lattice site i, Xi and pi denote the 
displacement and momentum of a harmonic oscillator at 
site i, and fii = '^^fii^ with fii^ — cl^c.^^. The third 
term, I^p, describes the coupling of dispersionless Ein- 
stein phonons to the electron occupation number n^. For 
doped cuprates or manganites, such a local interaction is 
expected to be a reasonable approximation as a result of 
screening. In the first term, the symbol (ij) denotes a 
summation over all nearest-neighbor hopping pairs 
and The parameters of the model are the hop- 

ping integral t, the phonon energy oj {h — 1), the el-ph 
coupling constant a, and the Coulomb repulsion [/ > 0. 
For J7 = 0, Eq. (^) is identical to the Holstein modelii^ 
As in previous work, we introduce the dimensionless cou- 
pling constant A — /{ujW), where W = 4iD is the bare 
bandwidth in D dimensions. We further define the pa- 
rameters uj = uj/t and U — U/t, and express all energies 
in units of t. Consequently, the independent parameters 
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of the model are aJ, A, and U. We shall see below that 
a very useful quantity is given by the polaron binding 
energy Ep = XW/2. Finally, throughout this paper, pe- 
riodic boundary conditions in real space are assumed. 

This work is exclusively concerned with the case of two 
electrons, neglecting the interaction between bipolarons 
which will definitely be present to some degree in real 
materials. Furthermore, we restrict our attention to two 
electrons with opposite spin, i.e., to the singlet bipolaron. 
A comparison of the singlet and triplet state has recently 
been given in Ref. 9i. 

A review of early work on the bipolaron problem 
can be found in the book of Alexandrov and MottiiS 
Here we focus the discussion on more recent devel- 
opments. The latter can roughly be divided into 
two classes depending on the methods employed: (a) 
variational approaches and (b) unbi- 
ased numerical studies using EDfi2*i2i2242ii2S variational 
diagonalization^SSiSiiSSiS^ the density-matrix renormal- 
ization group (DMRG)^ and QMC32a Except for the 
QMC study of de Raedt and Lagendijk,"?-- all work was re- 
stricted to the ground state. Moreover, even their QMC 
results were reported only for a single, low tempera- 
ture. This motivates our study of temperature effects 
in Sec. |IV| 

While ED and DMRG studies were obtained on clus- 
ters with twOfiSiSiiS^ four,^^ six,^'' eight^°i^^ or twelve 
sites 1 ^ the variational methods of Bonca et alm^ and of 
Refs. Inll2ll.'ll4llfll(fil are only weakly influenced by 
finite-size effects. An important disadvantage of ED and 
DMRG is the fact that the phonon Hilbert space has 
to be truncated, so that these methods can not easily 
be used to study the adiabatic (lU <C 1) and/or strong- 
coupling (A 3> 1) regime. In contrast, no such limitations 
are imposed on QMC and most variational methods. 

Although de Raedt and Lagendijk only consid- 
ered the adiabatic limit ZJ ~ 0, similar to other 
authors, ■'iiSiiiii^* their method can also be applied for 
finite phonon frequency.^* Moreover, it may be gener- 
alized to include dispersive phonons. Recently, an ex- 
tended Holstein model with long-range el-ph interac- 
tion has been investigated by Bonca and Trugman.'^" De 
Raedt and Lagendijk also considered long-range Coulomb 
interaction, while most other authors only took into ac- 
count the local Hubbard- type interaction given in Eq. 
except for Zhang et ali^ who have omitted this term 
in their DMRG calculations. Finally, we would like to 
point out that bipolaron formation in a model with Jahn- 
Teller modes — as present in the perovskite manganites — 
has been studied by Shawish et alm^ 



III. TRANSFORMED HAMILTONIANS 

The basis of both the variational approach and the 
QMC method presented below is the unitary transfor- 
mation H = vHv'^ of the Hamiltonian with v = 



exp(i^-j ^ijhiPj) (see Ref. 0). The result is 

H = -t^ 4cj.,e'^'(^'''^^')P'+P (2) 
^ « ■' 

K 

+ ^ hjX.,{Lj^ij - ady ) -I- ^ v,jfiihj - ^ TT-i , 

ij ij i 

^v^^^^— ^^^^^^ 

^cp ^ee 

with 

I 

As discussed in Ref. '7, henceforth also referred to as I, the 
extended transformation v takes into account nonlocal 
lattice displacements, which are essential for a correct 
description in the regime w < 1 . 

Similar to I, for the QMC method, we resort to the 
standard Lang-Firsov (LF) transformation"^^ with vq = 
exp(i7 fiiPi). Here 7 = ^ XW/uj has been chosen such 
that the el-ph coupling term I^p in Eq. cancels. The 
transformed Hamiltonian then takes the form 

Ho = -t^cLc^.,e'^(^^-w)+P 

^^^^^^^^^^^^^ 

A'o 

+ {U -2Ep)J2f^^^f^^l-'^Ep. (4) 



/ 

Hence, in contrast to the polaron problem,-i the el-el in- 
teraction term, resulting from the canonical transforma- 
tion, does not vanish but instead combines with the Hub- 
bard term. 



IV. QUANTUM MONTE CARLO 

The derivation of the QMC algorithm for the bipo- 
laron problem is very similar to the one-electron case^ 
and we shall therefore focus on the differences occurring. 
Moreover, we restrict the discussion to one dimension. 

A. Partition function 

We set out to calculate the partition function Z = 
g-/3ffo^ with ^0 given by Eq. To this end, we first 
notice that the last term in Hamiltonian Q is a constant 
and can therefore be neglected during the QMC simula- 
tion. Using the standard Suzuki- Trotter decomposition 
we obtain^ 
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where 13 = (fceT)"^ and At = (i/L. Inserting L complete 
sets of phonon momentum eigenstates and splitting up 
the trace into a bosonic and a fermionic part we findi 

Zl = Tr t j dpidp2 ■ ■ ■ dpL {pi\U\p2) ■■■ {pl I ^ bi> , 

(6) 

where dpr = Ilidpi.r: and lim^^oo -2^ = Z.^^ Since 
the phonon contribution to U is identical to the single- 
electron problem/ we can again integrate out the coordi- 
nates X. Upon defining Dp ~ dpidp2 ■ ■ - dpL the partition 
function becomes 



Zl / Vpwh Wf, 



(7) 



with C = [27r/(wAr)]^'^, 



(8) 



Here ifo.r is obtained from Kq [Eq. by replacing pi 
(pj) with pi^T- [pj^r)- The bosonic action has the form 



N 



(9) 



i=l 



with Pi — {pi^i, ■ ■ -Pi^l) and a tridiagonal L x L matrix 
A defined by 



A 



1 



Ll 



2 wAr2 



A 



1 



2wAt2 



(10) 



As pointed out in I, the representation of S'b given in 
Eq. permits us to introduce the so-called principal 
component representation discussed below. 

To evaluate the fermionic trace we choose the two- 
electron basis states 

[\l)^\i,j)^c\^c]^\Q) , i,j-l,...,iv} , (11) 

where we have introduced a combined index I running 
from 1 to N"^ in one dimension. We begin with the con- 
tribution of the kinetic term [Eq. |0J] . It follows that 
the tight-binding hopping matrix, denoted as k, has di- 
mension N"^ X N"^ . The exponential of the transformed 
hopping term can be written as^ 



(12) 



where 



We would like to emphasize that the random variables p 
merely enter the diagonal matrix Z?, while the N"^ x N"^ 
matrices Vr and k are fixed throughout the entire MC 
simulation. Thus, in total, we have 



n = WDrKD\.Vr ; 



(15) 



and the fermionic trace is calculated as 

which is identical to the sum over the diagonal elements 
of the matrix Q, in the basis . 

B. Observables 

The first observable of interest is the kinetic energy of 
the electrons defined as 

= -< E <4S.> = -2t5:(c,V^.,e"^(P--^^^-)) , (17) 

where we have exploited spin symmetry. Following the 
same steps as in the derivation of the partition function 
we get 



(gt^S^.,.) =2^'y I?p«;be'^(P-^-f-^)Trf(17c|^c^.^). (18) 
Writing out explicitly the fermionic trace we obtain 

(19) 



i'j' 



and the kinetic energy finally becomes 

= -2tZZ' I Vpw^Y.T."^''^"'^'"'''''^ {J,j'\^\hf) 



iv) j' 



(20) 

In addition to i?k, we shall also consider the correlation 
function 

p(<5) = E(n,Tn,+5^) , <5 = 0,l,...,iV/2-l. (21) 



A simple calculation leads to 



{Dr),v =dn'{dn,,,idn,,,i + Sn,^,iSn,,,i)e'^^P^^^+P^^-'> (13) p{5) ^ Z j u;b ^ (i, i (5| |z, z + <5) . (22) 



is diagonal in the basis lll|l . 

The second contribution to the matrix Q, in Eq. © 
comes from the effective el-el interaction term / [Eq. 
in terms of the diagonal matrix 



(14) 



Finally, we would like to point out that other observables, 
such as the total energy and the momentum distribution 
(cj,^Cj.^), may also be measured within the current ap- 
proach, while correlation functions such as {fiiXj) or the 
quasiparticle weight cannot be determined accurately^ 
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C. Principal components and reweighting 

We make use of the principal component representa- 
tion and the reweighting procedure, which have been dis- 
cussed in detail in I. Defining the principal components 

~ A^^^pi, in terms of which Sh [Eq. ©] takes a Gaus- 
sian form which can be sampled exactly/ allows to per- 
form calculations that are free of any autocorrelations 
between successive phonon configurations. In combina- 
tion with the reweighting, every new phonon configura- 
tion is accepted, and measurements can be made after 
each sweep through the N x L space-time lattice. The 
reweighting refers to the use of the purely bosonic weight 
Wb in the QMC simulation, while all the influence of 
the electrons and their interaction with the phonons — 
contained in W{ — is treated exactly as part of the observ- 
ables. 



D. Numerical details and performance 

The most significant difference between the present cal- 
culations and the one-electron case in I is the dimension 
of the matrices involved. While for one electron all ma- 
trices have size N x N — N being the extension of the 
ID lattice under consideration — here the dimension is 
X N"^. Clearly, this restricts calculations with re- 
spect to the number of lattice sites, especially in higher 
dimensions D > 1 where N'^ i-^ N'^^ . The total nu- 
merical effort for the current approach is proportional to 
N^^L. In contrast, the one-electron algorithm^ displays 
the same dependence cx N^^L as the determinant QMC 
method of Blankenbecler et a\m^ for the many-electron 
case, which can be reduced to N'^^L by employing the 
checkerboard breakup of the hopping matrix«^ The in- 
crease in required computer time for the bipolaron results 
from the fixed number of electrons. Recently, a grand- 
canonical version of the one-electron algorithm, also with 
a computer time ~ N'^^L, has been applied to study the 
dependence of polaron formation on carrier density in the 
spinless Holstein model,'^^ For the bipolaron problem, we 
shall see below that the present algorithm allows one to 
study lattices of reasonable size N < 14, for a wide range 
of the parameters oJ, A and U. In particular, we can 
obtain accurate results in the adiabatic regime ZU < 1. 

Let us briefly compare our method to other QMC ap- 
proaches to the HH bipolaron. The method of de Raedt 
and Lagendijk^^ is based on an analytic integration over 
the phonon degrees of freedom, leading to a model with 
retarded el-el interaction. Similar to our approach, it 
employs a Suzuki- Trotter approximation and gives re- 
sults at finite temperatures. For simplicity, de Raedt 
and Lagendijk only considered the adiabatic limit ZU = 0, 
in which there are no retardation effects. The numerical 
effort grows as L^, but is virtually independent of the 
system size, so that simulations can be carried out even 
for large clusters in three dimensions. However, it is not 
clear how a small but finite phonon frequency ZU < 1 will 



affect the computer time. 

Macridin et al,^^ used the diagrammatic QMC method 
to study two electrons on a 25 x 25 lattice. Although 
their approach does not rely on the Suzuki- Trotter de- 
composition, it is limited to zero temperature, and sta- 
tistical errors increase noticeably for ZU < 1. Moreover, 
the accuracy also decreases for large values of A and/or 
U, whereas we shall see in Sec. I VII that we can easily 
study the strong el-ph coupling regime also for J7 > 0. 

In I, we announced the possibility of reducing the nu- 
merical effort for the present method by exploiting the 
translational invariance of the model. To this end, the 
basis states ((TT|l would have to be replaced by states 
{\k, A)} with total quasimomentum k, and with the two 
electrons separated by a distance A. A similar idea has 
been used by KornilovitchSi^ for a single electron. In 
one dimension, the use of the basis {|fc, A)} would reduce 
the size of the matrices in the algorithm from N'^ x N'^ to 
N X N. However, in the course of the simulation, we had 
to evaluate the matrix product over r [Eq. ifTCjl ] for each 
allowed value of k. In total, we could therefore reduce 
the numerical effort by a factor N. The major drawback 
of using the reduced basis in momentum space is that 
it significantly complicates the program code. Conse- 
quently, in this work, we have restricted ourselves to the 
straight-forward extension of the one-electron algorithm 
presented in I. 

Finally, the minus-sign problem, which has been men- 
tioned in I, also exists here. However, as for one elec- 
tron, it quickly diminishes with increasing system size, 
and does therefore not conceivably affect simulations. 

V. VARIATIONAL APPROACH 

Although the method can easily be applied also in 
higher dimensions, we wish to keep the notation sim- 
ple and therefore restrict the derivation to D = 1. The 
approximation consists of the use of a zero-phonon basis 
after the extended unitary transformation, which leads to 
^op — [Eq. ||2Jl]. Furthermore, neglecting the ground- 
state energy of the oscillators, we also have P = 0, so 
that 

H = K + L, (23) 
with the transformed hopping term 

K = -ieff E = E ^(^) cLc,. (24) 

{ij)cr ka 

and e{k) = —2 tcS^ka^'^^i^)- Here, the effective hop- 
ping amplitude is given by^ 

ioff = -^e-3^'(^'-^~^')', (25) 

^ 5 

where 6 = ±1 in one dimension, z is the number of 
nearest neighbors, and rotational invariance has been ex- 
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ploited. For two electrons of opposite spin, the interac- 
tion term ||2J) simplifies to 

/oc = 2tio - [/ + 2 ^ v^jhi^fiji (26) 

if we use Vij = v\j-i\ and fiia-'hja = for i 7^ j. The 
two-electron eigenstates of the Hamiltonian H23|l have the 
form 

m = J2dp4-picU\o) ■ (27) 

Here we have suppressed the phonon component which 
is simply given by the ground state of N free harmonic 



oscillators. The states (|27ll may be written as 

\^'^) = ^T.^'''''T.'^A4+n\o) , (28) 

where the Fourier transform 

d = Fd (29) 

with Fip = e^^'P/^/N has been employed. The normal- 
ization of Eq. (|?7jl reads (V'fcl i^k) = J2p MpP- 

The expectation value of the transformed hopping term 
with respect to the states defined by Eq. (|27|l becomes 
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= J2\dpme{p)+s{k-p)] = -iUsd^Tkd. 
p 

In the last step we introduced vector notation, defined Tk — F diag[cos(p) -I- cos(fc — p)] /2 and used Eq. (|29|l . The 
expectation value of the interaction term is best computed in the real-space representation (|28|l . We find 

I ij j'j" W ' 



= {2vo - C/) E I'^'l' + § E"j+'.^-|^'l' = (2z;o - U)Sd + 2dWd 

3,1 



where the diagonal matrix Vij — SijVi has been intro- 
duced. The minimization with respect to d yields the 
eigenvalue problem 

(-4teff Tk + 2V)d^ [Eo ~2vo + U)d. (30) 

The vector of coefficients d and thereby the ground state 
are determined by minimizing the ground-state energy 
Eq through variation of the displacement fields 7^ . In the 
present work, we use the unconstrained nonlinear opti- 
mization routine fminsearch from the MATLAB package, 
together with several different starting points, including 
the simple LF result and random values of the 7^ . This 
ensures reproducible results even for a large number of 
variational parameters. 

In contrast to the local LF transformation, this proce- 
dure takes into account displacements of the oscillators 
not only at the same but also at the sites surrounding the 
two electrons. This represents a physically much better 
ansatz to describe the extended state which exists for 
weak el-ph coupling and/or strong Coulomb repulsion. 
Similar to the one-electron problem, we shall refer to the 



result obtained from the above variational method by re- 
placing 7ij with jSij as the Holstein Lang-Firsov (HLF) 
approximation. 

VI. RESULTS 

Before we turn to the results, we would like to review 
briefly the physics of the one-dimensional HH bipolaron 
as it emerges from existing work (see Sec. EJ. In the 
absence of Coulomb repulsion, the two electrons form a 
bound state for any A > 0. A crossover from an extended 
state, also called a large bipolaron, to a small bipolaron — 
with both electrons occupying the same site — is observed 
at a critical coupling strength Ac. The value of Ac is de- 
termined by the competition between the different terms 
in the Hamiltonian Similar to the one-electron case, 
for small phonon frequencies, the crossover takes place 
when the gain in potential energy due to bipolaron for- 
mation overcomes the loss in kinetic energy. While the 
former can be estimated in the atomic limit as AEp (see, 
e.g., Ref. 0), the latter is given here by = it — the ki- 
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TABLE I: Conditions for the existence of different singlet 
bipolaron states (see text). Here "wc" and "sc" denote weak 
coupling and strong coupling, respectively. 



[7 = 



Large bipolaron 


Small bipolaron 


A < 0.5 


A > 0.5 


or 


and 


2sjEvl^ < 1 


2^Ep/uj > 1 



U>0 


Two 


Intersite 


Small 


polarons 


bipolaron 


bipolaron 


U > 2Ep (wc) 
U > AEp (sc) 


U < 2Ep (wc) 
U < iEp (sc) 


[/<2Sp 



netic energy of the two electrons at A = 0. Since A can 
also be written in the form A = ^{AEp/W), we expect 
Ac = 0.5. For larger phonon frequencies cU 1, the lat- 
tice energy plays an important role, and gives rise to the 
additional criterion > 1 for the existence of a 

small bipolaroni^ 

For U > 0, a state with two weakly bound polarons is 
stable for weak enough el-ph interaction. Interestingly, 
starting from a small bipolaron a cross over to an intersite 
bipolaron — with the two electrons being localized most 
likely at neighboring lattice sites — takes place at a critical 
value {7cfi^*^^^ This state has been shown to have a 
much smaller effective mass than an on-site bipolaronjS^ 
and may therefore exist as a mobile quasiparticle in real 
systems. Phase diagrams of the intersite bipolaron have 
been reported in one^iS^ and two dimensions. For uJ = 
1, the region where such a state exists in one dimension 
is quite accurately described by J7 < 2Ep at weak el-ph 
coupling, and hy U < iEp at strong coupling!^ 

If 2Ep ^ U the two electrons can overcome the on-site 
repulsion and form a small bipolaron. As we shall see 
below, the competition between Hubbard repulsion and 
attractive interaction due to the electron-lattice coupling 
depends critically on the phonon frequency. A summary 
of the different bipolaron states together with (approx- 
imate) conditions for their existence is given in Tab.lil 

We would like to mention that very interesting physics 
can also be deduced from the bipolaron band dispersion 
(see Ref. i and references therein). Although the lat- 
ter can be studied in principle within our variational ap- 
proach, the subtle effects originating from the retarded 
nature of the el-ph interaction^ could not be addressed 
in a satisfactory way. 



A. Quantum Monte Carlo 

To eliminate the error ~ (Ar)^ due to the Suzuki- 
Trotter approximation, we extrapolate the QMC results 
to At = 0. In contrast, this error is expected to be 
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FIG. 1: (a) Normalized kinetic energy Sk [Eq. IjHl^ ] from 
QMC as a function of the el-ph coupling A for different values 
of the phonon frequency ZJ and the Hubbard repulsion U. 
(b) Dependence of Et^ on the cluster size N. Here and in 
subsequent figures QMC results have been extrapolated to 
At = (see text), errorbars are suppressed if smaller than 
the symbol size, and lines are guides to the eye. 



relatively large (on the order of a few percent) in the 
calculations of de Raedt and Lagendijk, due to the use of 
a rather small number of Trotter slices (L = 32 at /3 = 5, 
so that At ^ 0.16; see Ref. l2§) . Here we have performed 
simulations for three different values of At, typically 0.1, 
0.075 and 0.05. The errorbars in the figures below are 
usually as small as the linewidth, and will not be shown 
if smaller than the symbols used. 

Owing to the increased numerical effort compared to 
the one-electron case^^ we shall only present results for 
N < 12. Fortunately, finite-size effects on, e.g., the ki- 
netic energy, are already very small for this cluster size, 
as illustrated by Fig. ^b) for the most critical parame- 
ters. As expected, the largest changes with N occur near 
the crossover to a small bipolaron.^ Similar behavior has 
been found for the correlation function p{S). 

We define the effective kinetic energy of the two elec- 
trons as 



= E^/{-At) . 



(31) 



In Fig. ^a) we depict i?k as a function of the el-ph cou- 
pling for different values of oJ and U, at /3t = 10. This 
value of the inverse temperature is twice as large as in 
previous work)2& yielding results sufficiently close to the 
ground state to reveal the effects of bipolaron formation. 

Figure^a) reveals a strong decrease of E^ near A = 0.5 
for small phonon frequencies and U = 0. With increasing 
ZU, the crossover becomes less pronounced, and shifts to 
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larger values of A. For the same value of ZU, the crossover 
to a small bipolaron is sharper than the small-polaron 
crossover in the Holstein model with one electron (see, 
e.g., I). The small but finite kinetic energy even for strong 
el-ph interaction is a result of undirected, internal mo- 
tion of the two electrons inside the phonon cloud. For 
a finite on-site repulsion U = A, remains fairly large 
up to A ~ 1 (for tU < 2.0), in agreement with the strong- 
coupling result Ac = 1 for C/ = 4 (see discussion in Ref.|3)- 
At even stronger coupling, the Hubbard repulsion is over- 
come, and a small bipolaron is formed. Again, we see that 
the critical coupling increases with increasing phonon fre- 
quency. We would like to mention that the kinetic energy 
has also been calculated by ED on clusters of up to twelve 
sites r^^i^^i^ but results for ZU < 1 were restricted in the 
accessible range of A. In the regime where ED is appli- 
cable, a very good agreement has been found with our 
QMC data. 

The nature of the bipolaron state is revealed by the cor- 
relation function p{5) defined in Eq. (|21|l , which gives the 
probability for the two electrons to be separated by a dis- 
tance (5 > 0, and therefore represents a direct measure for 
the size of the bipolaron. Clearly, we have the sum rule 
"^sPi^) = 1- As pointed out, e.g., by MarsigHoji^ the 
phonon frequency determines the degree of retardation 
of the el-ph interaction, and thereby sets the maximal 
allowed distance between the two electrons compatible 
with a bound state. In the sequel, we shall focus on the 
most interesting case of small phonon frequencies, which 
has often been avoided in previous work for reasons out- 
lined in Sec. Ull 

Figure l^a) shows p{5) as a function of A for U = Q. 
Starting from the noninteracting state (A = 0) with p = 
1/iV, we see a pronounced increase of p(0) near A = 0.5. 
For large A > 2, we have p(0) ~ 1 and p{5) « for 5 > 0, 
characteristic for the on-site bipolaron. The decrease of 
the spatial extent of the bipolaron with increasing el-ph 
interaction is better illustrated in the inset of Fig. I^Ja), 
where we depict p as a function of 5. For finite on-site 
repulsion [/ = 4, an extended bipolaron state is stabilized 
for small A [Fig. |2Ib)], while a small bipolaron is found 
for A = 2. Additionally, we see that for A = 1, the 
electrons are most likely to occupy neighboring lattice 
sites [intersite bipolaron, see also inset in Fig.|2Ib)]. 

As pointed out earlier, a crossover from a small to an 
intersite bipolaron to two weakly bound polarons takes 
place as a function of the Hubbard interaction. Since 
the latter competes with the retarded el-ph interaction, 
the phonon frequency is expected to be an important pa- 
rameter. In Fig. |31 we show the kinetic energy and the 
correlation function p{5) as a function of U . We have 
fixed the el-ph coupling to A = 1. Starting from a small 
bipolaron for {7 = [see Fig.|21[a)], the kinetic energy in- 
creases with increasing Hubbard repulsion, equivalent to 
a reduction of the effective bipolaron massi^^^S^ Although 
the crossover is slightly washed out by the finite temper- 
ature in our simulations, there is a well-conceivable in- 
crease in i5k up to J7 w 4, above which the kinetic energy 




FIG. 2: Correlation function p(5) [Eq. I^TJ] from QMC as 
a function of el-ph coupling A for different values of 5. [(a) 
(7 = 0, (b) U = A\. Inset: Correlation function p{5) as a 
function of the distance 5 of the electrons for different A. 



begins to decrease again. The increase of originates 
from the breakup of the small bipolaron, as indicated 
by the decrease of p{Q) in Fig. Efb). Close to ?7 = 4, 
the curves for p(0) and p{l) cross, and it becomes more 
favorable for the two electrons to reside on neighboring 
sites. The intersite bipolaron only exists below a critical 
Hubbard repulsion Uc- As discussed at the beginning of 
this section, the latter is given by Uc = 2i?p (i.e., here 
Uc = 4) at weak el-ph coupling, and by Uc = 4:Ep at 
strong coupling. For an intermediate value A = 1 as 
in Fig. O the crossover from the intersite state to two 
weakly bound polarons is expected to occur somewhere 
in between, but is difficult to identify exactly from the 
QMC results. 

Figure 121 further illustrates that the crossover becomes 
steeper with decreasing phonon frequency. In the adia- 
batic limit tU = 0, it has been shown to be a first-order 
phase transition, whereas for ZJ > retardation effects 
suppress any nonanalytic behavior. At the same C/, E\, 
increases with uj since for a fixed A, the bipolaron be- 
comes more weakly bound. For the same reason, the 
crossover to an intersite bipolaron — showing up in Fig. |31 
as a crossing of p(0) and p(l) — shifts to smaller values of 
U. 

Let us now consider the effect of temperature. While 
the kinetic energy shows a similar dependence as in the 
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FIG. 3: (a) Normalized kinetic energy and (b) correlation 
functions p(0), p(l) from QMC as a function of the Hubbard 
repulsion U for different values of the phonon frequency ZU. 



one-electron case — with the crossover being smeared out 
at high temperatures — it is much more interesting to 
look at In Figs. QJa) - (c) we plot p(S) at differ- 

ent temperatures, for parameters corresponding to the 
three regimes of a large, small and intersite bipolaron, 
respectively. 

Large bipolaron. For the parameters chosen (U — 0, 
A ~ 0.25), the two electrons are most likely to occupy 
the same site, but the bipolaron extends over a distance 
of several lattice constants [Fig. ^a)]. Clearly, in this 
regime, the cluster size TV = 12 used here is not com- 
pletely satisfactory, but still provides a fairly accurate de- 
scription as can be deduced from calculations for iV = 14 
(not shown). Nevertheless, on such a small cluster, no 
clear distinction between an extended bipolaron and two 
weakly bound polarons can be made. As the tempera- 
ture increases from f3t — 10 to pt = 1, the probability 
distribution broadens noticeably, so that it becomes more 
likely for the two electrons to be further apart. In partic- 
ular, for the highest temperature shown, p(0) has reduced 
by about 30 % compared to f3t = 10. 

Small bipolaron. A different behavior is found for the 
small bipolaron, which exist at stronger el-ph coupling 
A = 1.0. Figure 0fb) reveals that p{6) peaks strongly at 
6 = 0, while it is very small for 5 > at low tempera- 
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FIG. 4: Correlation function p{S) from QMC as a function of 
5 for different inverse temperatures P, N = 12 and ZJ = 0.4. 



tures. Increasing the temperature, we observe that p{S) 
remains virtually unchanged up to pt = 3. Only at very 
high temperatures there occurs a noticeable transfer of 
probability from 6 = to 6 > 0. At the highest tem- 
perature shown, (3t = 0.5, the two electrons have a non- 
negligible probability for traveling a finite distance 5 > 
apart, although most of the probability is still contained 
in the peak located at (5 = 0. 

Intersite bipolaron. Finally, we consider in Fig. 01c) 
the intersite bipolaron, which has been found above for 
U = 4 and A = 1.0 [Fig.|2tb)]. At low temperatures, p{S) 
takes on a maximum for 5 = 1. For smaller values of (3t, 
the latter diminishes, until at /3t = 1, the distribution is 
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completely flat, so that all 6 are equally likely. 

The dilTerent sensitivity of the bipolaron states to 
changes in temperature found above can be explained 
by their different binding energies. The latter is given 
by AEo = - 2E''q \ where e'^^^ and e'^^^ denote the 
ground-state energy of the model with one and two elec- 
trons, respectively. These quantities can be calculated 
using the present method as well that presented in I. 

Generally, the thermal dissociation is expected to oc- 
cur at a temperature such that the thermal energy k^T = 
{/3T)~^ becomes comparable to AEq, in accordance with 
our numerical data. The large and the intersite bipolaron 
are relatively weakly bound as a result of the rather small 
effective interaction Ucs m U — 2Ep (see discussion in 
Ref.H). The binding energies are AEq « -(0.32 ± 0.08)t 
and —(0.28 ± 0.08)t, respectively, so that we expect a 
critical inverse temperature (31 « 2.5-5?^ in agreement 
with Figs. EJa) and (c). In contrast, the small bipo- 
laron in Fig.^b) has a significantly larger binding energy 
AE « —(3.43 ± 0.09)i, and therefore remains stable up 
to pt « 0.3. 

Since the thermal dissociation of intersite bipolarons 
has been proposed to explain the activated dc conduc- 
tivity in the paramagnetic state of the manganites;^ we 
would like to comment on the relation of our findings 
to this theory. Instead of the Holstein-type model used 
here, Alexandrov and BratkovskySi^ argue in favor of a 
model with long-range el-ph interaction, and assume that 
charge carriers are O p holes rather than Mn d holes, so 
that the double-exchange mechanism does not come into 
playi^ An intersite bipolaron in their theory is formed 
by two holes residing on neighboring oxygen ions. Fur- 
thermore, they also include a nearest-neighbor Coulomb 
repulsion V ~ Ep between electrons. In the present case, 
the latter would, most importantly, reduce the binding 
energy of the intersite state, thereby leading to a lower 
temperature for dissociation. For sufficiently large V, in- 
tersite bipolaron formation is expected to be completely 
suppressed. A closer investigation of this issue in the 
framework of the Holstein-Hubbard model may be car- 
ried out using a generalization of the present method. 

In total, a quantitative comparison of our numerical 
results to the work of Refs. and to the 3D mangan- 
ites, appears to be not justified due to the simplifications 
made and the different model studied in the present work. 



B. Variational approach 

While the above QMC approach is limited to finite 
temperatures and relatively small clusters, the varia- 
tional method of Sec. yields ground-state results on 
much larger systems. It becomes exact in several limits. 
First, for A = (i.e., no cl-ph coupling), we obtain the 
exact solution 7^- — for all i, j. Second, as cxd, 
no phonons can be excited so that the use of a zero- 
phonon basis is justified. Similarly, in the classical limit 
uj = 0, the phonons do not have any dynamics, and the 



0.6 



IW 



0.4 - 



0.2 



\ 


1 ' 1 


1 

N = 25,U = " 


1 Q ^■ 
T ^\ ^ 


N 

\ 

\ 


— m = 0.4 




■y. \ 
' T\ \ 


03 = 1.0 


\ 




(0 = 2.0 
03 = 4.0 
<^ 03 = 0.4, HLF 




\ 


^ CO = 4.0, HLF 




\ 

\ 

\ 

\ 





2 

X 



FIG. 5: Variational results for the normalized kinetic energy 
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FIG. 6; Variational results for the normalized kinetic energy 
-Bk and the correlation functions p(0), p(l) as a function of 
the on-site repulsion U. 



variational determination of the displacement fields al- 
lows one to obtain exact results for any A. In contrast, 
the HLF approximation (see Sec. IVj generally overesti- 
mates the displacement of the lattice at a given site in 
the presence of an electron, even for ut — 00. Finally, the 
variational approach becomes exact in the nonadiabatic 
strong-coupling limit A,w — *■ 00. Since the two-electron 
problem is diagonalized exactly without phonons, the 
above statements hold for any value of the Hubbard re- 
pulsion U. 

To scrutinize the quality of the variational method, we 
started by comparing the ground-state energy for [7 = 
as a function of el-ph coupling for different values of UJ, 
to the most accurate approach currently available in one 
dimension, namely the variational diagonalization^ We 
find a good agreement over the whole range of A. As 
expected from the nature of the approximation, slight 
deviations occur for cJ < 1. 
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Despite the success in calculating the total energy — 
being the quantity that is optimized — one has to be care- 
ful not to overestimate the validity of any variational 
method. To reveal the shortcomings of the current ap- 
proach, we show in Fig. [S] the normalized kinetic en- 
ergy Ek — tcff [see Eqs. and as a function 
of el-ph coupling, and for different U. We have chosen 
iV = 25 to ensure negligible finite-size effects. In princi- 
ple. Fig. [S] displays a behavior similar to the QMC data in 
Fig.^a). There is a strong reduction of i?k near A = 0.5 
for U = 0.4, which becomes washed out and moves to 
larger A with increasing phonon frequency. Compared 
to the exact QMC results in Fig.^a), the crossover to a 
small bipolaron is much too steep in the adiabatic regime, 
regardless of the fact that the variational results are for 
T = 0. This is a common defect of variational methods. 
Moreover, for U — 0.4-2.0, the variational kinetic energy 
is too small above the bipolaron crossover compared to 
the QMC data, while for oJ ~ A, the decay of with 
increasing A is too slow. 

The reason for the failure is the absence of retarda- 
tion effects, which play a dominant role in the forma- 
tion of bipolaron states. The increased importance of 
the phonon dynamics — not included in the variational 
method — for the two-electron problem leads to a less 
good agreement with exact results than in the one- 
electron casei In particular, our variational results over- 
estimate the position of the crossover (Fig. |SJ| compared 
to the value Ac = 0.5 expected in the adiabatic regime. 
Nevertheless, the method represents a significant im- 
provement over the simple HLF approximation, due to 
the variational determination of the parameters 7ij . This 
is illustrated in Fig. [S] where we also show the HLF re- 
sult i?k — exp{—Ep/Lu) for oJ = 0.4 and 4.0. In con- 
trast to the variational approach, the HLF approxima- 
tion yields an exponentially decaying kinetic energy for 
all values of the phonon frequency. Whereas such behav- 
ior actually occurs in the nonadiabatic limit w — > oo, the 
situation is different for small uj [see Figs, ^^a) andE]. 
Thus the HLF approach cannot reproduce the physics of 
bipolaron formation for small and intermediate phonon 
frequencies, while the variational method presented here 
accounts qualitatively for the dependence on the phonon 
frequency. 

Next, we wish to study the influence of Coulomb repul- 
sion U. Similar to Fig.|31 we take A = 1, so that an on-site 
bipolaron state is formed at J7 = 0. For small phonon 
frequency ZJ = 0.4, Fig. reveals a sharp crossover near 
U = 2.5, i.e., at a smaller value of U than in the QMC 
results of Fig. O the reason again being the neglect of 
the retarded nature of the effective el-el interaction. As 
in the QMC results, the Coulomb repulsion breaks up 
the on-site bipolaron, leading to an increase of the ki- 
netic energy. Moreover, the curve for p{l) peaks at the 
crossover point, indicating the existence of an intersite 
bipolaron in this regime. A similar picture is found for 
larger phonon frequency Zj = 4, also shown in Fig. El 
although the changes with increasing U are much more 
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FIG. 7: Upper panel: variational results for the effective 
interaction Ueft{S) (see text) and the correlation function p{S) 
(inset) as a function of the el-el distance S. Lower panel: 
variational lattice distortions 7^ as a function of 5. 



gradual than for uj — 0.4. 

Finally, we report in Fig. [71 (upper panel) the effective 
interaction Ucs{S) between the two electrons as a func- 
tion of their relative distance S, given by vs [Eq. P])]. 
We have chosen uj — 0.4 and U — 4, the same param- 
eters as in Fig. |2Ib). For A = 0.75, the finite Coulomb 
repulsion stabilizes two weakly bound polarons, as illus- 
trated by the results for p{d) shown in the inset of Fig. [7\ 
While UcS is repulsive (positive) for (S = 0, the two elec- 
trons can form a bound state by traveling a finite dis- 
tance 1 < (5 < 4 apart. This is still true for A = 1, for 
which the HLF approach yields J7eff(0) = U — 2Ep — 0. 
Nevertheless, the two electrons experience an attractive 
interaction and form an extended bipolaron. Finally, for 
even stronger coupling A = 1.25, the phonon-mediated el- 
el interaction has overcome the on-site repulsion, so that 
Ucs{S = 0) < 0. At the same time, the size of the bipo- 
laron has collapsed to a single site. This crossover is also 
well visible in the lower panel, which displays the vari- 
ationally determined lattice distortions 75. It is worth 
mentioning that the values of Ucs{0) in Fig. [3 are larger 
than the strong-coupling prediction U—2Ep for all values 
of A considered. This may be attributed to the overesti- 
mated bipolaron binding energy in the atomic limit. 

As pointed out in several places, the shortcomings of 
the variational approach presented here are a result of the 
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missing dynamical phonon effects. The present approach 
may be further improved by making an ansatz for the 
eigenstates of the untransformed Hamiltonian of the 
form 

ij p 

x(j(%t{^W} + 4^).t{^(^)}) 1^,,), 

with defined as in Eq. Hll(l . two canonical transfor- 
mations depending on the displacement fields 7,^^'' and 
7^^^ (see Sec. lIII|) . and additional variational parameters 

Sp \ Thereby, one can take into account lattice dis- 
tortions not centered at the sites of the electrons, which 
become important as tU — > 0, and which reproduce to 
some degree the effect of retardation. 



Our results underline the importance of the phonon dy- 
namics, which has often been neglected in previous work. 
Moreover, we have presented for exact results for the ef- 
fect of temperature on the bipolaron state in the im- 
portant adiabatic regime. Thermal dissociation of bipo- 
larons has been observed at temperatures where the ther- 
mal energy becomes comparable to the binding energy. 

Two interesting open issues are the effect of nearest- 
neighbor Coulomb interaction, as well as that of dimen- 
sionality. While the latter cannot easily be addressed 
with the current approach, one may instead extend the 
promising work of Ref. to finite phonon frequencies. 

Finally, we have proposed a variational ansatz based 
on a canonical transformation with variational parame- 
ters. The latter represents a significant improvement over 
standard approximations. In particular, it qualitatively 
accounts for the dependence on the phonon frequency. 



VII. CONCLUSIONS 



We have studied the Holstein-Hubbard bipolaron with 
quantum phonons by extending a quantum Monte Carlo 
method previously developed for the Holstein polaroni 
In its present form, the method is limited to one- 
dimensional clusters. However, in contrast to other ap- 
proaches, it allows to perform accurate calculations also 
for small phonon frequencies and finite temperatures. 

We have studied the dependence of bipolaron forma- 
tion on the phonon frequency and the Hubbard repulsion. 
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